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FINAL  REPORT 

AFOSR/DARPA  Contract  F49620-91-C-0059 
Mathematical  Library  Software 
for  Applications  of  Parallel  Supercomputers 


1.  OVERVIEW 


In  this  work  we  have  developed  application  software  for  a  variety  of 
advanced  scientific  problems  on  a  variety  of  parallel  supercomputer  platforms. 
The  supercomputer  platforms  used  in  this  work  include  the  IBM  PVS,  IBM  SPl 
and  IBM  SP2,  Cray  T3D,  and  Silicon  Graphics  Power  Challenge.  The  IBM  PVS, 
located  in  our  research  facility  at  Princeton  University,  has  been  used  for  many  of 
the  research  studies  reported  herein.  It  has  32  processors  and  the  demonstration 
of  good  scaling  behavior  on  this  system  has  been  one  of  our  primary 
achievements. 

Applications  of  our  work  have  ranged  from  micro-lithography  to  chemical 
kinetics  to  fluid  turbulence.  Among  the  significant  breakthroughs  accomplished 
with  support  of  this  contract  include  the  first  aerial  image  of  a  full-scale  4Mb 
DRAM  chip  (in  fact,  the  IBM  4Mb  memory  chip)  which  initially  required  about 
one  day  on  our  IBM  PVS  to  calculate  but,  with  algorithmic  improvements,  can 
now  be  done  in  less  than  one  quarter  of  that  time.  We  have  also  made  significant 
breakthroughs  in  the  development  of  parallel  lattice  gas  methods,  especially 
lattice  BGK  methods,  which  enable  natural-grid  solution  of  complex  problems. 
In  fluid  turbulence,  breakthroughs  have  included  the  first  direct  simulation  of  the 
Kolmogorov  k"5/3  inertial  range  spectrum  and  the  development  of  hyperviscous 
simulation  techniques  that  allow  parallel  calculation  at  Reynolds  numbers  that 
rival  the  highest  Reynolds  numbers  observed  in  nature,  namely  geophysical 
flows. 
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A  variety  of  parallel  application  software  has  been  developed  and  has 
been  communicated  to  a  broad  variety  of  collaborators  in  industry,  universities, 
and  government  laboratories. 

In  the  remainder  of  this  report,  we  give  a  list  of  the  papers  published  with 
support  of  this  contract  and  a  brief  summary  of  the  results  obtained  in  each  of 
these  papers. 
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Number  of  Students  Participating  in  the  AFOSR/DARPA  and  the  Number  of 
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A.  S.  Ghandi 
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Y.  Zhang 


Ph.D4/95 
Ph.D  expected  6/96 
Ph.D  expected  8/96 
Ph.D  expected  8/95 
Ph.D  5/95 
Ph.D  expected  8/95 
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II.  BIBLIOGRAPHY 


Papers  Acknowledging  Support  of 
Contract  F49620-91-C-0059 


Interaction  of  Surface  Waves  with  Turbulence:  Direct  Numerical  Simulations 
of  Turbulent  Open-Channel  Flow  (Vadim  Borue,  Steven  A.  Orszag  and  Ilya 
Staroselsky,  J.  Fluid  Mech.  (1995),  vol.  286,  pp.  1-23. 

We  report  direct  numerical  simulations  of  incompressible  unsteady  open- 
channel  flow.  Two  mechanisms  of  turbulence  production  are  considered:  shear 
at  the  bottom  and  externally  imposed  stress  at  the  free  surface.  We  concentrate 
upon  the  effects  of  mutual  interaction  of  small-amplitude  gravity  waves  with  in- 
depth  turbulence  and  statistical  properties  of  the  near-free-surface  region. 
Extensions  of  our  approach  can  be  used  to  study  turbulent  mixing  in  the  upper 
ocean  and  wind-sea  interaction,  and  to  provide  diagnostics  of  bulk  turbulence. 

Numerical  Simulation  of  Weakly  Compressible  Kolmogorov  Flow  with 
Kinetic  Model  (Y.  H.  Qian  and  S.  A.  Orszag),  Conference  Proceedings  & 
Lecture  Notes  in  Physics,  Vol.  1,  "Second  International  Conference  on 
Computational  Physics,"  ICCP-2,  Sept.  13-17,  1993,  Beijing,  China, 
International  Press  (1995). 

Kolmogorov  flow  has  been  used  to  study  the  transition  to  turbulence.  The 
motoviation  of  this  paper  is  to  investigate  this  flow  with  a  simple  kinetic  scheme. 
We  observe  numerically  the  -5/3  law  of  inverse  cascade  of  energy  spectrum  at 
low  wavenumbers  and  a  steeper  than  -3  law  at  high  wavenumbers  when  the  flow 
is  forced  at  small  scales. 

Lattice  BGK  Models  for  Thermohydrodynamics  (Y.  H.  Qian  and  S.  A.  Orszag), 
Conference  Proceedings  &  Lecture  Notes  in  Physics,  Vol.  1,  "Second 
International  Conference  on  Computational  Physics,"  ICCP-2,  Sept.  13-17, 1993, 
Beijing,  China,  International  Press  (1995). 

Lattice  BGK  models  are  proposed  for  thermohydrodynamics.  The 
theoretical  values  of  the  sound  speed,  the  shear  viscosity  and  the  conductivity 
are  verified  by  numerical  simulations.  One-dimensional  shock-tube  problems  is 
used  to  test  the  models. 
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Pattern  Formation  in  Phase  Transition  (Y.  H.  Qian  and  S.  A.  Orszag), 
Conference  Proceedings  &  Lecture  Notes  in  Physics,  Vol.  1,  "Second 
International  Conference  on  Computational  Physics,"  ICCP-2,  Sept.  13-17, 1993, 
Beijing,  China,  International  Press  (1995). 

The  dynamics  of  liquid-gas  phase  transition  is  studied  in  this  paper  by 
using  simple  kinetic  models.  The  analytical  results  are  verified  numerically  in 
one,  two  and  three  dimensions. 

Parallel  Computation  for  Fluid  Dynamics  (Y.H.  Qian  and  S.  A.  Orszag), 
Conference  Proceedings  &  Lecture  Notes  in  Physics,  Vol.  1,  "Second 
International  Conference  on  Computational  Physics,"  ICCP-2,  Sept.  13-17, 1993, 
Beijing,  China,  International  Press  (1995). 

There  is  a  growing  interest  in  parallel  algorithms  for  computational  fluid 
dynamics.  In  this  paper,  we  describe  some  simple  kinetic  models  which  are 
capable  of  reproducing  the  hydrodynamic  behaviours  with  a  high  efficiency. 

Scalings  in  Diffusion-Driven  Reaction  A  +  B  C:  Numerical  Simulation  by 
Lattice  BGK  Models  (Y.H.  Qian  and  S.  A.  Orszag),  Conference  Proceedings  & 
Lecture  Notes  in  Physics,  Vol.  1,  "Second  International  Conference  on 
Computational  Physics,"  ICCP-2,  Sept.  13-17,  1993,  Beijing,  China, 
International  Press  (1995). 

We  are  interested  in  one  application  of  lattice  BGK  models  to  the 
diffusion-driven  reactive  system  A  +  B  C,  which  was  investigated  by  Galfi  and 
Racz  with  an  asymptotic  analysis  and  by  Chopard  and  Droz  with  a  cellular 
automaton  model.  The  lattice  BGK  model  is  free  from  noise  and  flexible  for 
applications.  We  derive  the  general  diffusion-reaction  equations  for  the  lattice 
BGK  models  under  the  assumption  of  local  diffusive  equilibrium.  The  4th  order 
terms  are  derived  and  verified  by  numerical  simulations.  The  purpose  of  this 
paper  is  to  compare  the  lattice  BGK  results  with  existing  results  before  we  apply 
the  models  to  more  complicated  systems.  The  scalings  concern  the  exponents  a , 
b  and  the  function  G  in  the  production  rate  of  C  component  R(x,f)  =  t~^G(xt 
"®).  We  find  the  same  values  for  a  =  1/6  and  b  =  2/3  as  Glfi  and  Racz  found  at 
long  time  limit.  A  Gaussian  function  for  G  is  numerically  obtained  for  the  first 
time.  Compared  with  the  asymptotic  analysis,  lattice  BGK  is  easy  to  apply  to 
cases  where  no  analytic  or  asymptotic  results  exist.  Compared  with  the  cellular 
automaton  model,  the  lattice  BGK  model  is  faster,  simpler  and  more  accurate. 
The  discrepancy  of  the  results  between  the  cellular  automaton  model  and  the 
lattice  BGK  models  for  the  exponents  come  from  the  role  of  the  intrinsic 
fluctuation.  Once  the  time  and  space  correlation  of  noise  is  given,  we  can 
incorporate  a  stochastic  term  in  our  models.  The  Schlogl  model  is  also  tested  in 
this  paper. 


4 


Recent  Advances  in  Lattice  Boltzmann  Computing  (S.  Succi,  Y.  H.  Qian  and  S. 
A.  Orszag),  submitted  for  publication,  1995. 

Recent  developments  in  the  theory  of  the  Lattice  Boltzmann  equation  are 
presented,  in  particular,  lattice  BGK  models  are  discussed  in  detail.  Various 
applications  are  summarized. 

Self-Similar  Decay  of  Three-Dimensional  Homogeneous  Turbulence  with 
Hyperviscosity  (V.  Borue  and  S.  A.  Orszag),  Phys.  Rev.  E,  51,  R856-R859  (1995). 

Numerical  simulations  of  the  Navier-Stokes  equations  with  h5q3erviscosity 
(_i[)/i+lA^(h=8)  show  that  periodic-box  turbulence  exhibits  self-similar  decay. 
The  inertial-range  energy  spectrum  has  the  scaling  law  Z(t)  2/3y^5/3  where  Z(t) 
is  the  energy  dissipation  rate  at  time  t .  The  total  energy  of  the  system  decreases 
as  1/ The  concept  of  constant  Reynolds  number  decay  is  introduced,  enabling 
us  to  perform  long  time  averages  and  reliably  measure  higher  order  correlation 
functions.  Comparisons  are  made  with  the  case  of  forced  turbulence  reported 
earlier. 

Numerical  Studies  of  Three-Dimensional  Kolmogorov  Flow  at  High  Reynolds 
Number  (V.  Borue  and  S.  A.  Orszag),  submitted  for  publication. 

High  resolution  numerical  simulations  (up  to  256^  modes)  are  performed 
for  the  three-dimensional  flow  in  a  periodic  box  of  size  L  =  2n  driven  by  the  large 
scale  force /y  =  FcosM  (Kolmogorov  flow).  High  Reynolds  number  is  attained  by 

solving  the  Navier-Stokes  equations  with  the  hyperviscosity  (-1)^+^A^  =  8).  It 

is  shown  that  the  mean  velocity  profile  of  the  Kolmogorov  flow  is  practically 

independent  of  Reynolds  number  and  has  the  "laminar"  form  aV  cos(x)  with 
nearly  constant  eddy  viscosity.  Nevertheless  the  flow  is  highly  turbulent:  the 
turbulent  kinetic  energy  is  nearly  equal  to  the  mean  flow  kinetic  energy.  The 
turbulent  intensities,  the  energy  dissipation  rate  and  different  terms  in  the  energy 
balance  equations  have  very  simple  structure  of  coordinate  dependences  a  + 
bcos(2x)  (here  a,  b  are  some  constants).  This  makes  Kolmogorov  flow  a  good 
testing  case  for  studying  turbulent  transport  models.  At  large  scales  the  flow  is 
anisotropic,  although  at  high  Reynolds  number  it  becomes  isotropic  at  small 
scales.  The  problem  of  local  isotropy  is  systematically  studied  by  measuring 
longitudinal  and  transverse  components  of  energy  spectra  and  cross-correlations 
of  different  components  of  velocity.  It  is  shown  that  pressure  plays  the  crucial 
role  in  making  flow  locally  isotropic.  It  is  demonstrated  that  nearly  a  decade  in 
wave-number  space  is  required  before  anisotropic  large  scale  flow  may  be 
considered  locally  isotropic. 
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Forced  Three-Dimensional  Homogeneous  Turbulence  with  Hyperviscosity  (V. 
Borue  and  S.  A.  Orszag),  Europhysics  Letters  29, 687-692  (1995). 

High  resolution  numerical  simulations  (up  to  256®  modes)  are  performed 
for  the  Navier-Stokes  equations  with  h5q)erviscosity  (-l)^+^A(h=8).  With  random 
white-in-time  forcing,  the  inertial-range  energy  spectrum  has  the  scaling  law 

1/k^  with  a  ~  1.85  ±  0.05  at  small  wave  numbers  k  with  a  flatter  region  near 
the  ultraviolet  cutoff.  The  transitional  region  is  characterized  by  a  sharp  change 
of  the  energy  transfer  mechanism.  Spectra  of  fluctuations  of  higher  order 
correlation  functions  are  measured.  At  inertial  range  scales  the  system  seems 
insensitive  to  the  power  h  of  the  hyperviscosity. 

Inverse  Energy  Cascade  in  Stationary  Two-Dimensional  Homogeneous 
Turbulence  (V.  Borue),  Phys.  Rev.  Lett.  72, 1475-1478  (1994). 

Direct  numerical  simulations  with  up  to  1024^  resolution  are  performed  to 
study  statistical  properties  of  the  inverse  energy  cascade  in  stationary 
homogeneous  two-dimensional  turbulence  driven  by  small  scale  Gaussian  white- 
in-time  noise.  The  energy  spectra  for  the  inverse  energy  cascade  deviate  strongly 
from  the  expected  kr^^^  law  and  are  close  (somewhat  flatter)  to  k~^.  The  reason 
for  the  deviation  is  traced  to  the  emergence  of  strong  vortices  distributed  over  all 
scales.  Statistical  properties  of  the  vortices  are  explored. 

Aerial  Image  of  3D  Phase-Shifted  Reticle  --  3D  FAIM  (E.  Barouch,  U. 
Hollerbach,  and  S.  A.  Orszag),  Proc.  SPIE  2197, 442, 1994. 

As  is  well  know,  the  mask  in  projection  printing  is  a  thin  (about  1  jim 
thick)  object,  composed  of  numerous  variable  domains.  The  various  domains 
contain  different  complex  refractive  indices.  These  refractive  indices  form  a 
discontinuous  "mask  function",  appearing  as  a  coefficient  in  the  Maxwell- 
Material  (MM)  equation  that  describes  the  transmission  of  the  electric  field 
through  the  3D  mask  domain.  The  right-hand  side  of  this  equation  involves  the 
gradient  of  the  dot-product  of  the  electric  field  with  the  gradient  of  the  log  of  the 
refractive  index.  Within  each  feature,  the  refractie  index  is  a  constant.  Therefore 
its  gradient  vanishes,  and  the  resulting  local  MM  equation  is  the  wave  equation. 
However,  at  a  feature  boundary  (chrome-quartz,  phase-1  —  phase-2,  etc.)  the 
complex  refractive  index  is  discontinuous,  resulting  in  a  highly  singular  delta- 
function-like  coefficient.  The  trailing  components  of  the  field  display  the 
standard  Maxwell  disoontinuity  conditions,  resulting  in  two  very  singular  terms, 
i.e.  a  derivative  of  a  delta  function  and  a  product  of  two  delta  functions.  Such  a 
source  term,  with  such  a  remarkable  oscillation,  will  be  the  origin  of  extreme 
instability  in  the  numerical  solution  of  the  equations,  no  matter  which  direct 
solution  algorithm  is  employed.  Accordingly,  the  derivation  of  a  preconditioner, 
namely  a  function  which  contains  most  of  the  unpleasantness  of  the  equation 
times  a  function  to  be  found,  is  absolutely  essential  to  obtain  accurate  solutions  in 
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a  reasonable  time.  Through  an  elaborate  uniformly  asymptotic  scheme  and  scale 
analysis,  we  have  derived  such  a  preconditioner.  The  solution  of  the  direct  MM 
equation  was  compared  with  the  solution  using  the  preconditioner,  and  it  gave 
us  the  confidence  that  we  were  on  the  right  track.  It  actually  means  that  one 
wishes  to  find  the  mask's  Green  function  which,  when  convolved  with  the 
electric  field  incident  on  the  mask,  yields  the  exiting  field  which  then  enters  the 
imaging  optical  system.  It  is  the  purpose  of  this  paper  to  demonstrate  the  results 
obtained  thus  far,  as  well  as  to  illustrate  its  implications  on  defocus  linewidth 
control. 

Automated  Determination  of  CAD  Layout  Failures  Through  Focus: 
Experiment  and  Simulation  (C.  Spence,  J.  Nistler,  E.  Barouch,  U.  Hollerbach, 
and  S.  A.  Orszag),  Proc.  SPIE  2197, 302, 1994. 

This  paper  describes  a  software  program  which  compares  the  simulated 
aerial  image  of  a  mask  pattern  with  the  original,  desired,  mask  pattern.  The 
program  uses  the  Fast  Aerial  Image  Model  (FAIM)  to  simulate  the  aerial  image. 
A  specific  intensity  contour  is  chosen  and  the  distance  from  the  contour  to  the 
desired  design  is  calculated.  Regions  where  the  distance  exceeds  a  specified 
tolerance  are  deemed  failures  and  flagged.  Corner-rounding  errors  are  handled 
differently  to  line  edge  position  errors.  The  threshold  intensity  used  can  be 
specified  by  the  user,  alternatively  a  'critical  feature'  may  be  defined  and  the 
threshold  set  to  the  intensity  value  required  to  print  it  on  size.  In  addition  to 
describing  the  program,  we  also  show  examples  of  how  the  aerial  image  contour 
predictions  compare  with  simulations  of  resist  profiles  and  actual  printed  resist 
images  for  the  case  of  an  SRAM  cell. 

Spectral  Exponents  of  Enstrophy  Cascade  in  Stationary  Two-Dimensional 
Homogeneous  Turbulence  (V.  Borue),  Phys.  Rev.  Lett.  71, 3967-3970,  December 
1993. 


Direct  numerical  simulations  with  up  to  4096^  resolution  are  performed  to 
address  the  question  of  universality  of  statistical  properties  of  the  enstrophy 
cascade  in  homogeneous  two-dimensional  turbulence  driven  by  large-scale 
Gaussian  white-in-time  noise.  Data  with  different  Reynolds  numbers  are 
compared  with  each  other.  The  energy  spectrum  is  foimd  to  be  very  close  to 
1  /fc3.  It  is  shown  that  the  primary  contribution  to  the  enstrophy  transfer  function 
comes  from  wave-number  triads  with  one  small  leg  and  two  long  ones, 
corresponding  to  wave  numbers  in  the  inertial  range. 
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Reynolds  Number  Dependence  of  Isotropic  Navier-Stokes  Turbulence  (Z-S 
She,  S.  Chen,  G.  Doolen,  R.  Kraichnan,  and  S.  A.  Orszag),  Phys.  Rev.  Lett.  70, 
3251-3254,  May  1993. 

Reynolds  number  dependence  of  turbulence  energy  spectra  and  higher- 
order  moments  of  velocity  differences  is  explored  by  numerical  integrations  of 
the  incompressive  Navier-Stokes  equation.  The  simulations  have  spatial 
resolutions  up  to  512^  and  cover  15  <  Rx  <  200,  where  Rx  is  the  Taylor 
microscale  Re5molds  number.  The  energy  spectra  collapse  when  scaled  by  the 
wave  number  kp  of  peak  dissipation  and  also  by  the  spectrum  level  at  kp  .  kp 

varies  with  Rx  in  accord  with  the  1941  Kolmogorov  theory.  High-order 
normalized  moments  of  velocity  differences  over  inertial-range  distances  exhibit 
an  R;\^-independent  variation  with  separation  distance.  Implications  of  these 
observations  are  discussed. 

Far-Dissipation  Range  of  Turbulence  (S.  Cheng,  G.  Doolen,  J.  Herring,  R. 
Kraichnan,  S.  A.  Orszag,  and  Z.  S.  She),  Phys.  Rev.  Lett.  70,  3051-3054,  May 
1993. 


Navier-Stokes  turbulence  at  low  Reynolds  number  {Rx  ~  15)  is  studied  by 
high-resolution  computer  simulation.  The  energy  spectrum  in  the  range  5kd  < 
lOkd  is  well  fitted  by  k  ^  exp{-ck/kd),  where  kd  is  the  Kolmogorov  dissipation 

wavenumber,  a  =  3.3  and  c=7.1.  High-order  spatial  derivatives  of  the  velocity 
field  exhibit  strong  intermittency,  associated  with  gentle  spatial  variation  of 
large-scale  structure  rather  than  with  sparse,  intense  small-scale  structures. 
Analysis  by  the  direct-interaction  approximation,  which  ignores  intermittency, 

gives  a  =  3,  c=ll. 

Vector  Aerial  Image  with  Off-Axis  Illumination  (E.  Barouch,  D.  C.  Cole,  U. 
Hollerbach  and  S.  A.  Orszag),  Proc.  of  the  SPIE 1927, 686, 1993. 

As  the  numerical  aperture  (NA)  in  optical  projection  systems  increases, 
the  vector  nature  of  the  projected  electric  and  magnetic  fields  becomes  more 
important.  Recent  advances  in  off-axis  illumination,  tilting  condenser  lenses,  and 
applying  spatial  filters,  as  well  as  phase-shift  masks,  hold  the  potential  of 
increased  depth  of  focus.  To  simulate  these  techniques  in  the  high-NA  regime, 
we  have  developed  new  fast  altorithms.  The  development  reported  here  builds 
on  our  recent  work  [1]  extending  scalar  aerial  imaging.  The  systematic  treatment 
of  [1]  has  now  been  applied  to  account  for  the  vector  nature  of  light.  The  optical 
projection  system  examined  here  involves  an  aplanatic  lens  system  with  high  NA 
and  finite  magnification  as  applied  to  both  i-line  and  deep-UV,  4x  and  5x, 
steppers.  Our  new  algorithms  owe  their  genesis  to  a  variety  of  subtle 
mathematical  properties  of  the  scalar  and  vector  aerial  image  formulae.  These 
new  algorithms  are  based  on  our  high  NA  aerial  image  model  juxtaposed  with 
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Yeung's  original  ideas  for  the  vector  field  image.  This  formulation,  which 
involves  non-uniform  grid  fast  Fourier  transform  techniques,  has  resulted  in  a 
tremendous  speedup  of  the  aerial  image  simulation  with  uniform  accuracy 
throughout  a  domain  of  lOOOp  x  1000|i  and  beyond,  with  the  potential  to 
simulate  an  entire  chip,  thus  providing  a  solution  to  the  "grand  challenge"  of  full- 
chip  design. 

Lattice  BGK  Models  for  the  Navier-Stokes  Equation:  Nonlinear  Deviation  in 
Compressible  Regimes  (Y.  Y.  Qian  and  S.  A.  Orszag),  Europhys.  Lett.,  21(3), 
255-259, 1993. 

In  this  letter,  a  nonlinear  deviation  from  the  Navier-Stokes  equation  is 
obtained  from  the  recently  proposed  LBGK  models,  which  are  designed  as  an 
alternative  to  lattice  gas  or  lattice  Boltzmann  equation.  The  classical  Chapman- 
Enskog  method  is  extended  to  derive  the  nonlinear-deviation  term  as  well  as  its 
coefficient.  Their  analytical  expression  is  derived  for  the  first  time,  thanks  to  the 
simplicity  of  the  LBGK  models.  A  numerical  simulation  of  a  shock  profile  is 
presented.  The  influence  of  the  correction  on  the  kinetics  of  compressible  flow  is 
discussed.  A  complete  analysis  of  the  thermodynamics  including  the 
temperature  will  be  presented  elsewhere. 

Direct  and  Large-Eddy  Simulation  of  Axisymmetric  Wakes  (A.  G. 
Tomboulides,  S.  A.  Orszag  and  G.  E.  Karniadakis),  AIAA  31st  Aerospace 
Sciences  Meeting,  Paper  93-0546,  January  1993. 

In  this  paper  we  present  direct  numerical  simulations  (Re  =25- 1,000)  and 
a  large  eddy  simulation  {Re  =  20,000)  for  flow  past  a  sphere,  which  is  a  protot5q3e 
axisymmetric  body.  A  new  spectral  element-Fourier  method  is  implemented  to 
discretize  the  incompressible  Navier-Stokes  in  cylindrical  coordinates,  and  a 
sub-grid  model  based  on  the  renormalization  group  method  (RNG)  is  used  in  the 
large  eddy  simulation.  We  identify  two  early  bifurcations  of  the  sphere  wake, 
the  first  at  Reynolds  number  212  leading  to  a  three-dimensional  time- 
independent  state,  and  the  second  at  Reynolds  number  between  250  and  285 
resulting  in  a  time-periodic  flow.  A  shear  layer  instability,  observed  in 
experiments,  is  also  resolved  accurately  at  Reynolds  number  1,000.  Our  results 
from  the  large  eddy  simulation  are  in  good  agreement  with  experimental  results 
both  in  terms  of  temporal  response  (frequency  spectrum)  as  well  as  in  terms  of 
transport  measures  (drag  coefficient). 

Modes,  Nodes,  and  Flow  Codes  (G.  E.  Karniadakis  and  S.  A.  Orszag),  Physics 
Today,  34-42,  March  1993. 

The  numerical  simulation  of  high  Reynolds  number  turbulent  flows  is  a 
grand  challenge  problem  that  will  tax  supercomputing  capability  for  decades  to 
come.  Massively  parallel  supercomputers  seem  the  best  hope  for  achieving 
progress  on  these  fundamental  problems  with  many  practical  applications. 
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Some  Basic  Challenges  for  Large  Eddy  Simulation  Research  (S.  A.  Orszag,  L 
Staroselsky,  and  V.  Yakhot),  in  "Large  Eddy  Simulation  of  Complex 
Engineering  and  Geophysical  Flows,"  Cambridge  Lfniv.  Press  (ed.  by  B. 
Galperin  and  S.  A.  Orszag),  pp.  55-78, 1993. 

This  paper  reviews  some  of  the  basic  features  of  LES  and  the  techniques 
used  to  derive  the  LES  equations.  We  discuss  the  conditions  for  self-consistency 
of  LES.  Then,  review  the  RNG  theory  of  the  Navier-Stokes  equations  and  the 
RNG-based  derivation  of  a  subgrid  model  for  LES.  We  discuss  the  effect  of  rigid 
walls  on  LES  and  illustrate  the  problematics  of  the  LES  modeling  of  complex 
flows  on  the  example  of  a  subgrid-scale  model  for  turbulent  combustion. 

On  the  Small-Scale  Dynamical  Behaviour  of  Lattice  BGK  and  Lattice 
Boltzmann  Schemes  (S.  Succi,  D.  d'Humieres,  Y.  H.  Qian,  and  S.  A.  Orszag),  J. 
Sci.  Comp.  8, 219-230, 1993. 

A  quantitative  assessment  of  the  validity  of  the  lattice  BGK  version  of 
Lattice  Boltzmann  schemes  is  presented  for  the  case  of  two-dimensional 
incompressible  turbulence. 

Eliminating  Spurious  Invariants  in  Lattice  Models  for  Hydrodynamics  (Y.  H. 
Qian  and  S.  A.  Orszag),  submitted  to  Complex  Systems,  to  be  published. 

Most  of  lattice  gas,  lattice  Boltzmann  and  lattice  BGK  models  suffer  from 
the  existence  of  so-called  spurious  invariants  which  are  usually  scalar  and 
transported  by  the  dynamics.  Their  influence  is  shown  in  the  pressure  term.  The 
presence  of  spurious  invariants  drives  the  solution  of  the  system  away  from  the 
Navier-Stokes  equation's  solutions.  Though  in  lattice  Boltzmann  or  lattice  BGK 
models,  the  initial  quantities  of  these  spurious  invariants  can  be  set  to  zero  in 
order  to  minimize  their  influence  on  dynamics,  the  production  of  such  invariants 
by  shock  wave  (strong  gradients)  or  boundary  conditions  is  not  yet  clear.  We 
think  that  the  existence  of  spurious  invariants  is  one  of  the  main  sources  of 
numerical  instabilities.  So  it  is  important  to  remove  the  spurious  invariants.  In 
this  paper,  we  will  be  able  to  eliminate  some  of  such  spurious  invariants  by 
replacing  the  propagation  step  by  a  smoothing  procedure.  We  rederive  the 
governing  equations  and  obtain  the  transport  coefficients.  Some  stagged 
invariants  can  be  completely  eliminated  if  the  smoothing  parameter  Q  is  set  to 
0.5.  Numerical  simulation  does  confirm  the  theoretical  results.  The  consequence 
of  removing  small-scale  spurious  invariants  is  meaningful  and  useful  for  subgrid 
modelling,  and  we  do  believe  that  the  smoothing  procedure  is  very  general  and 
that  unphysical  behaviours  can  be  reduced. 
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